H2 formation via non-Born-Oppenheimer hydrogen migration in photoionized ethane

Neutral H2 formation via intramolecular hydrogen migration in hydrocarbon molecules plays a vital role in many chemical and biological processes. Here, employing cold target recoil ion momentum spectroscopy (COLTRIMS) and pump-probe technique, we find that the non-adiabatic coupling between the ground and excited ionic states of ethane through conical intersection leads to a significantly high yield of neutral H2 fragment. Based on the analysis of fingerprints that are sensitive to orbital symmetry and electronic state energies in the photoelectron momentum distributions, we tag the initial electronic population of both the ground and excited ionic states and determine the branching ratios of H2 formation channel from those two states. Incorporating theoretical simulation, we established the timescale of the H2 formation to be ~1300 fs. We provide a comprehensive characterization of H2 formation in ionic states of ethane mediated by conical intersection and reveals the significance of non-adiabatic coupling dynamics in the intramolecular hydrogen migration.

repetition rate of 1 kHz. The laser beam was focused onto the molecular beam using a mirror with a focal length of 75 mm and its polarization was oriented vertically with respect to the time-of-flight (TOF) direction. The process of molecular ionization and fragmentation induced by the laser pulses was investigated by the cold target recoil ion momentum spectroscopy (COLTRIMS). The produced electrons and ions were detected in coincidence by two timeand position-sensitive detectors in a uniform magnetic field and electric field configuration.

Supplementary note 2: Surface hopping molecular dynamics simulation and frequency analysis
In order to validate the interpretation of the non-adiabatic effect in the experiment, we performed fewest-switches surface hopping (FSSH) molecular dynamics (MD) simulation on the CASSCF(5,5)//6-31G * level to determine the branching ratio of the dissociative H 2 formation channel. The MD simulation includes non-adiabatic couplings for the lowest three doublet cationic states in the adiabatic representation, denoted as D 0 , D 1 and D 2 , which are correlated with the cationic 2 E g and 2 A 1g states in the D 3d symmetric geometry. The D 0 state is correlated with the diabatic 2 A 1g state in the dissociation limit, and the energy barrier for H 2 formation reaction initiated from Franck-Condon geometry C 2 H + 6 is 0.68 eV. As shown in Figure 4 of the main text, the FSSH simulation results in a 67% branching ratio of the dissociative H 2 formation channel in the 2 E g state, and the calculated branching ratio of the bound parent ion channel in the 2 A 1g state is 31%, showing consistency with the experimental results. The FSSH molecular dynamics simulation results indicate that the non-adiabatic coupling between 2 E g and 2 A 1g states is crucial for the H 2 formation dynamics of ethane cation C 2 H + 6 . Furthermore, we performed frequency analysis for the MD trajectories, and investigate its relationship with the frequency distribution for the CI geometries with H 2 formation character. The result, together with the frequency analysis of the experimentally measured C 2 H 2+ 4 + H 2 yields oscillation, are shown in Figure 5 in the main text. Continuous wavelet transformation (CWT) method are used to obtain the time and frequency distribution I CWT (ω, t) from velocity autocorrelation function where i is index of atoms, and ⟨· · · ⟩ represents the average over molecular dynamics trajectories. The CI geometry shown in Fig. 1c in the main text is optimized by the quantum chemistry packages Molpro [1] initialized from the transition state geometry of the H 2 formation reaction C 2 H + 6 → C 2 H + 4 + H 2 . We define the characteristic frequency distribution of CI geometry to be where ω i is the i-th vibrational mode frequency of the optimized geometry around the conical intersection, and the vibrational modes are shown in Supplementary Figure 2. The peak time for the overlap between I CWT (ω, t) and I CI (ω) corresponds to the characteristic time when C 2 H + 6 reaches the CI geometry. Supplementary note 3: Fitting procedure details for above-threshold ionization spectrum The above-threshold ionization (ATI) spectroscopy contains typical comb structures, which are spaced by the photon energy of the laser pulses. These combs arise from electrons tunneling from the two highest occupied molecular orbitals (HOMO and HOMO-1) and their energy offsets reflect the differences in ionization potentials between these orbitals. Consequently, the comb structures observed in the ATI spectra can serve as a signature to disentangle contributions from multiple orbitals in ATI spectroscopy. In the low-intensity region, polyatomic molecules, including ethylene and other molecules mentioned in the literature [2], exhibit high contrast ATI combs. For ethane, the nonadiabatic coupling between the ionic ground and excited states plays a key role in the formation of H 2 , as depicted in Figure 1 of the main text. The ATI combs in coincidence with the parent ion and fragment consist of two components corresponding to the cationic ground state 2 E g (resulting from HOMO ionization) and the excited state 2 A 1g (resulting from HOMO-1 ionization), respectively. Considering the present peak laser intensity, direct tunneling ionization (TI) from HOMO-2 can be neglected, as the ionization rate for HOMO-2 is two orders of magnitude smaller than that for HOMO, according to molecular Ammosov-Delone-Krainov (ADK) calculations. To disentangle the contributions from the 2 E g and 2 A 1g electronic states, we applied a global fitting procedure to reconstruct the observed ATI spectra from different channels, as shown in Supplementary Figure 3. To minimize the influence of the Freeman resonance (FR), the fitting procedure ignores the spectrum below 1 eV. Additionally, extending the energy range beyond 7 eV does not affect the fitting results, as the ATI combs are completely smeared out. Thus, the fitting is performed within the energy range of 1∼7 eV in the ATI spectrum. In panels (a) and (b), the laser intensity is 88 TW/cm 2 , while in panels (c) and (d), the laser intensity is 82 TW/cm 2 . All the spectra are fitted with two series of Gaussian functions (see Supplementary . The solid cyan line represents the fitted curve, which reproduced the measurement (black) very well. The fitted components, denoted as Comp. 1 and Comp. 2, correspond to the contributions from HOMO and HOMO-1 ionization, respectively. The calculated ATI spectras obtained using the Coulomb-corrected strong-field approximation (CCSFA) approach are displayed in the green and red shadows, which corresponds to the ionization from HOMO and HOMO-1 orbitals of ethane, respectively.
The fitting function we used is as follows: where Y j are the measured ATI spectra signals in coincidence with the parent ion (j = 1) and fragment ion (j = 2) of ethane, which are normalized with respect to the total electron counts. B j % is the ratio of the 2 E g state in the  To estimate the influences of background signals on the fitting procedure, the events from false coincidence and noise in the raw ATI spectra need to be considered. In our measurement, the count rates for electrons and ions are approximately 0.1/pulse and 0.075/pulse, respectively, with a total detection efficiency of around 50%. Therefore, the estimated false coincidence rate is 3%. The average noise level in each measurement is below 5%. To determine the total error of the fitting results, we added a random noise level of 10% to the total counts of each energy bin in the raw ATI spectrum and generated 1000 sampled spectra. The same fitting procedure was then applied to these spectras, and the error of the fitting results was obtained and shown in Supplementary Table  1 for two laser intensities. We utilized ethylene as a benchmark for the analysis of ATI spectra. Ethylene exhibits an energy gap of approximately 2 eV between its first excited and ground ionic states. In the fitting process, we considered two components and observed two overlapping energy combs from the fitting, as shown in Supplementary  Figure 4, suggesting that ionization from the HOMO is dominant in ethylene.

Supplementary note 4: Reconstruction of PMDs corresponding to 2 E g and 2 A 1g states
The measured and calculated PMDs of ethylene are shown in Supplementary Figure 5c and d, respectively, demonstrating good agreement and proving the validity of CCSFA calculation. Figure 5e and f show the measured raw PMDs in coincidence with the parent ion and H 2 formation channels for ethane, while the calculated PMDs for the HOMO ( 2 E g state) and HOMO-1 ( 2 A 1g state) ionization are presented in Figure 5h and j, respectively. A discrepancy between the PMDs from the raw measurement and the calculations can be observed. The distinct angular distributions cannot be easily interpreted using the simple explanation commonly applied to diatomic molecules [3][4][5]. Based on our molecular dynamics (MD) simulations and fitting approach, we suggest that both the 2 E g and 2 A 1g states The measured PMDs in coincidence with the parent ion and H 2 formation channels respectively. g, h The reconstructed PMDs corresponding to the 2 E g and 2 A 1g states of ethane respectively. i, j The calculations compared to the PMDs in (g) and (h). The laser intensity is 76 TW/cm 2 for ethylene while 88 TW/cm 2 for ethane. In (g) and (h), the fanning out stripes are located within the pink semicircle for the 2 E g state. In (i) and (j), the arm-like structure is dominant within the pink semicircle. The jet-like structures in the ATI spectra are emphasized by the two solid blue curves.
contribute to the PMDs of the parent ion and H 2 formation channels due to non-adiabatic coupling dynamics through conical intersections (CIs). According to the fitted values of B 1 % and B 2 %, the ATI spectrum of ionization from the HOMO and HOMO-1 can be reconstructed, as discussed in Supplementary note 3. It shows the measured PMDs of the parent ion and H 2 formation channels are linear combinations of PMDs corresponding to 2 E g (HOMO) and 2 A 1g (HOMO-1) states with different relative ratios, which can be described as following: Thus the PMD corresponding to HOMO and HOMO-1 can be reconstructed from the measured PMDs once we know B j %. However, directly applying the same method used for ATI to reconstruct the PMD solely from the HOMO or HOMO-1 is not feasible, as the angular-dependent TI yields have different relative contributions of the two components for various electron emission angles. To overcome this challenge, we divided the PMD into six angular regions and obtained the relative yields of the two components by fitting each integrated ATI spectra (formalized as Supplementary  Equation 6). The value of six is chosen simply because of the limit of the statistics of electron counts. Since the PMDs of the parent ion and H 2 formation channels both contain contributions from the 2 E g and 2 A 1g states with different relative ratios, the PMD exclusively from the 2 E g state can be reconstructed by subtracting the normalized PMD of H 2 formation channel from that of parent ion channel, where the normalization ratio is determined by B 1 %.  show the angular distributions of ATI peak (12-photon channel) for 2 E g and 2 A 1g states. The numbers of the maxima are labeled next to the position of the peaks. All the curves have been symmetrized.
The reconstructed PMDs of the 2 E g and 2 A 1g states, based on the fitted values (Supplementary Table 2), are shown in Supplementary Figure 5g and i. This method can also be applied to the case of 82 TW/cm 2 laser intensity, enabling the reconstruction of the PMDs shown in Fig. 3c and d in the main text. These reconstructed PMDs demonstrate good agreement with the calculations based on the CCSFA approach. The characteristic features, such as the jet-like distributions of the ATI ring (within the solid blue curves) and the low-energy structures (within the pink semicircle), are well reproduced by the calculations, validating the accuracy of the reconstructed PMDs (in Supplementary Figure  6). In the case of the 2 E g state, the fanning out structure contains of six stripes (11-photon channel), and the jet-like structure within ATI ring (12-photon channel) shows seven nodes, which are consistent with the calculations. For the 2 A 1g state, the jet-like structure within the ATI (12-photon channel) contains seven nodes, matching the calculation. At a higher laser intensity of 88 TW/cm 2 , the arm-like structure becomes less pronounced compared to 82 TW/cm 2 , mainly due to the electron energy from the 11-photon channel in ATI being close to 0 eV. The quantitative agreement between the extracted and simulated PMDs strongly confirms the reliability of the fitting procedure.
The transition amplitude M p in strong-field approximation (SFA) in the form of saddle point approximation can be written in the atomic units as where A and E are the vector potential and electric field of laser, |ψ 0 ⟩ is the initial bound state of the electron, |p+A(t)⟩ is the Volkov state, S(t) = t −∞ dt ′ [p+A(t ′ )] 2 /2 + I p t is the action. Saddle point solutions t s can be obtained by solving the saddle point equation The above expressions are known as the SFA, which states that once the electron is ejected, Coulomb interaction between the electron and the residual ion is ignored. The CCSFA theory adds corrections to the SFA by considering the effect of Coulomb field to the trajectories and action of ionized electrons [6] Here p is the asymptotic momentum disturbed by the Coulomb field, S ′ (t s ) is the corrected action which reads The first integral in the right hand of Supplementary Equation 10 denotes the sub-barrier action, where p is the undisturbed asymptotic momentum, the upper limit t r is the real part of t s . The second integral denotes the action obtained after the electron leaving the tunnel exit. The initial state of molecules |ψ 0 (t)⟩, within the fixed-nuclei approximation, can be written as a linear combination of atomic orbitals (LCAO). Here we consider the symmetry properties of molecular orbitals under the mirror reflection perpendicular to C-C bond. The C 2 H 4 molecule has D 2h symmetry and the C 2 H 6 molecule has D 3d symmetry, and their HOMO and HOMO-1 orbital symmetry properties are shown in Supplementary Figure 5(a,b) and Figure 1, respectively. The resulting LCAO orbitals can be written as combinations of atomic orbital pairs. In each pair, the two centers of the atomic orbitals are symmetric about the origin just as in a homonuclear diatomic molecule, such that pairs of atomic orbitals form symmetric (γ = 1) and antisymmetric (γ = −1) combinations: where R a denotes the relative nuclear coordinates, the subscript a denotes different atom pairs and γ can be 1 or −1 depending on the symmetry. We show the contributions from the different centers of the molecules to the molecular orbitals in Supplementary Table 3.
After substituting the wavefunction Supplementary Equation 11 into Equation 9 according to the dressed modified molecular SFA, the transition amplitudes can be rewritten as  where f γ,a (p, R a ) = 2 cos (p · R a /2) 2i sin (p · R a /2) is the interference term between two atomic centers in one atom pair. For some molecules, not only the HOMO but also the lower-lying orbitals can substantially contribute to the ionization process. Hence, we include both the HOMO and HOMO-1 orbitals in our calculations. In Supplementary  Table 4, we show the ionization potentials of different orbitals. The difference in ionization potential between the HOMO and HOMO-1 orbitals for ethylene is 2 eV, indicating a dominant contribution of HOMO to the ionization. Moreover, our calculations reveal that the transition amplitude of HOMO-1 is at least two orders of magnitude smaller than that of HOMO for ethylene. On the other hand, for ethane, the difference in ionization potential between the HOMO and HOMO-1 orbitals is only 0.7 eV. In this case, the contributions from both the HOMO-1 and HOMO orbitals become comparable.
In Supplementary Figure 5e to h, the prominent orbital-dependent features are observed in the low-energy region within the pink semicircle for the 2 A 1g state, specifically the fanning out structure for the 2 E g state and the armlike structure for the 2 A 1g state. These orbital-dependent structures can be understood based on the parity of the molecular orbitals' wave functions. Within the fixed-nuclei approximation, the numerically calculated molecular initial state for ethane can be given as a linear combination of atomic orbitals (LCAO) for one carbon-carbon (C-C) pair and three hydrogen-hydrogen (H-H) pairs. As shown in Supplementary Figure 7, if we separate the contributions of each atomic pair, we can see that the arm-like structure is prominent in the PMDs of the H-H pairs but absent in the PMD of the C-C pair. The parity of the total wave function, and hence the corresponding PMD, depends on the combination of different atomic pairs. For the 2 E g state, only two pairs of H-H contribute, and they are combined in odd parity (see Supplementary Table 3). As a result, the contributions of the two H-H pairs to the ionization process cancel each other out, and the total PMD is mainly determined by the contribution of the C-C pair. This is confirmed by the calculated PMD from the C-C pair (Supplementary Figure 7a and e), which reproduces the reconstructed PMD of the 2 E g state well. However, for the 2 A 1g state, the three pairs of hydrogen atoms are combined in even parity. Therefore, the total PMD is determined by both the C-C pair and the H-H pairs, resulting in the appearance of the arm-like structure in the very low-energy region.
The strong-field ionization depends on the alignment of molecules. Actually, the simulated energy spectra and PMDs presented in the manuscript are the averaged results among different orientations to mimic the experimental conditions in which the molecules are presumably randomly aligned. In Supplementary Figures 8 and 9, the simulated ionization yields and PMDs for the 2 E g and 2 A 1g states of ethane are shown for the various angles between the C-C axis and the laser polarization. According to the calculations, the ionization yield of the 2 E g state exhibits a strong dependence on the molecular orientation. It reaches a maximum at an angle of 45 degrees between the C-C axis and the laser polarization, while minima are observed at 0 and 90 degrees. Therefore, the orientation-averaged results, as presented in the manuscript, are predominantly determined by the ionization of molecules oriented around 45 degrees. In contrast, the ionization yield of the 2 A 1g state shows a much weaker dependence on the orientation angle. The different behaviors of the two stated root in the different constitutions of their wave functions. The ionization of the 2 E g state is mainly influenced by the C-C atomic pair, making it sensitive to molecular orientation. But for the 2 A 1g state, both C-C pair and H-H pairs contribute to the ionization. The three H-H pairs orienting in different directions smear out the orientation effect. In addition, the pattern of the PMDs for both states exhibits a weak dependence on molecular orientation. The critical features used to distinguish the two states, i.e., the fanning out structure for the 2 E g state and the armlike structure for the 2 A 1g state, are prominent at nearly all orientation angles. This is because the pattern of the PMD at low energy is primarily determined by the symmetry of the molecular orbital, which is retained at different orientations. Based on these calculations, it is expected that the molecular alignment has a minor effect on the reconstructed results of the two states under the considered experimental conditions.  We established that the non-adiabatically coupled 2 E g and 2 A 1g cation states can both contribute to the parent ion and H 2 formation channels. In Supplementary Table 5, we provided the method to extract the branching ratios and listed the important numbers of electron events for two laser intensities. Moreover, we qualitatively show that the sub-cycle TI yield ratios of the HOMO-1 orbital (N A /(N P +N D ) has an enhancement with increasing the laser intensity from 82 TW/cm 2 to 88 TW/cm 2 . These values are shown in the last row of Supplementary Table 5. The error bars of all branching ratios from TI in Figure 4 of the main text are obtained by propagating the error of B 1 % and B 2 % shown in Supplementary Table 1.
We obtained the branching ratios from the FR ionization. As discussed in the main text, the peak e appears for both channels, and the peak a is only pronounced for the dissociative H 2 formation channel. The branching ratios of dissociation channel for the two electronic states can be also derived by the ratio of counts between two channels without the fitting procedure used in Figure 2 of the main text. The numbers of the electron events from the integration are shown in Supplemetary Table 6. The error range of branching ratios can be obtained by propagating the error of the counts of peak e or peak a. By comparing the results induced by TI and FR, we can claim that the obtained branching ratios are solely determined by the non-adiabatic coupling dynamics between the ground and excited cationic states. 88 TW/cm 2 peak e (Branching ratio) peak a (Branching ratio) Parent NeP=307 (NeP/Ne = 20(∓5)%) NaP=29 (NaP/Na = 11(±10)%) Dissociation NeD=1246 (NeD/Ne = 80(±5)%) NaD=234 (NaD/Na = 89(∓10)%) Total Ne=NeP + NeD=1553 Na=NeP + NaD=263 82 TW/cm 2 peak e (Branching ratio) peak a (Branching ratio) Parent NeP=5564 (NeP/Ne = 28(∓4)%) NaP=764 (NaP/Na = 19(±8)%) Dissociation NeD=14337 (NeD/Ne = 72(±4)%) NaD=3241 (NaD/Na = 81(∓8)%) Total Ne=NeP + NeD=19901 Na=NeP + NaD=4005 decreases from 1.0 eV to 0.5 eV, and at the same time (Supplementary Figure 10c), the yield of this band clearly increases for time delays larger than 1000 fs. These features of the band can reflect the neutral H 2 formation dynamics. The appearance time of this band (approximately 700 fs) corresponds to the moment when the wave-packet starts to arrive at the CI and the H 2 formation channel becomes active, consistent with the results presented in the main text for channel (1). Due to the large-scale motion of the hydrogen atom during the non-adiabatic dynamics along the cation, the KER of this band is smaller compared to the direct Coulomb explosion channel from the dication. As the delay time increases, the distance between C 2 H + 4 and H 2 groups becomes larger, resulting in an even smaller KER (Supplementary Figure 10c). The neutral H 2 formation dynamics underlying the time-dependent features of the (C 2 H + 4 + H + 2 ) channel are overall consistent with channel (1). However, the yield of channel (2) is much smaller than channel (1), so the main focus in the main text is on channel (1). At the early time before the wave-packet arrives to the CI, we observed another channel C 2 H + 6 → CH + 3 + CH 3 (3). This channel shows opposite time-dependent dynamics compared to channel (1) (C 2 H 2+ 4 + H 2 ), as shown in Supplementary Figure 11. Channel (3) originates from the dissociation of higher excited states of the cation. It exhibits a clear enhancement starting from 500 fs and starts to decrease after 1300 fs. Furthermore, the increasing counts of channel (1) and the decreasing counts of channel (3) are at the same level, suggesting that the excitation of C 2 H + 6 is an alternative pathway before the wave-packet arrives at the CI. After passing through the CI, the ionization of C 2 H + 4 becomes more significant relative to channel (3). This result indicates that the initial wave-packet, after ionization, can be promoted to higher excited states induced by the probe pulse before the C 2 H + 4 is reformed.